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Abstract In this paper, a fine numerical algorithm, namely the "clean numerical 
simulation" (CNS), is proposed to accurately simulate the propagation of micro- level 
uncertainty of chaotic dynamic systems. The chaotic Hamiltonian Henon-Heiles sys- 
tem for motion of stars orbiting in a plane about the galactic center is used as an 
example to show the validity of the CNS. It is found that, due to the sensitive depen- 
dence on initial condition (SDIC) of chaos, the micro-level inherent uncertainty of the 
position and velocity of a star transfers into the macroscopic randomness of motion. 
Thus, chaos is a bridge from the micro-level inherent uncertainty to the macroscopic 
randomness in nature. This provides us a new explanation to the SDIC of chaos from 
the physical viewpoint. 

1 Introduction 

Using high performance digit computers, a lots of complicated problems in science, 
finance and engineering have been solved with satisfied accuracy. However, there 
exist some problems which are still rather difficult to solve even by means of the most 
advanced computers. One of them is the propagation of micro-level uncertainty of 
chaotic dynamical systems. 

It is well-known that all numerical simulations are not "clean" : there exist more 
or less numerical noises such as truncation error and round-off error, depending on 
numerical algorithms. In most cases, such kind of numerical noises are much larger 
than the micro-level uncertainty of dynamic systems under consideration, so that the 
micro-level uncertainty is completely lost in the numerical noise. This becomes more 
serious for chaotic dynamic systems, which have the sensitive dependence on initial 
conditions (SDIC), i.e. small difference of initial condition leads to great difference 
of numerical simulations of chaotic solution at large enough time. Thus, very fine 
numerical algorithms need be proposed so as to accurately simulate the propagation 
of micro-level uncertainty of chaotic dynamic systems. This is the motivation of this 
article. 

Physically, many chaotic dynamic systems contain micro-level uncertainty in ini- 
tial conditions. For example, let us consider the famous chaotic dynamic system 
governed by Lorenz equation, which is based on Navier-Stokes equation under the 
so-called continuum- assumption of fluid. In physics, any fluids are made of atoms 
or molecules, and thus not a continuum. So, the concepts in fluid mechanics such 
as velocity, density, pressure, as so on, are defined in the statistic meanings: strictly 
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speaking, all of them have statistic fluctuations and thus are inherently uncertain. 
Although these micro-level uncertainty might be much less than numerical noises of 
traditional numerical algorithms, they may propagate together with the chaotic so- 
lution and, due to the SDIC of chaos, have a great influence on chaotic solution at 
long enough time. Such kind of propagation of micro-level uncertainty has important 
physical meanings, and may deepen our understandings about the SDIC of chaotic 
dynamic system from physical and statistic view-points. 

In this article, a kind of fine numerical algorithm, called "clean" numerical simula- 
tion (CNS), is proposed to accurately simulate propagation of micro- level uncertainty 
of chaotic dynamic systems. Here, the word "clean" means that the truncation error 
and round-off error are much less than the micro-level uncertainty so that the nu- 
merical noises can be neglected in a long enough interval of time for the propagation 
of uncertainty. A chaotic Hamiltonian system proposed by Henon and Heiles [7J is 
used to show its validity. The basic ideas of the so-called clean numerical simulation 
(CNS) are given in §2, followed by the investigation of the micro-level uncertainty of 
the system in §3 and its propagation in §4 from statistical points of view. Conclusions 
and discussions are given in §5. 



2 The numerical algorithm of the CNS 

Henon and Heiles [7J proposed a Hamiltonian system of equations 

x = —x — 2xy, (1) 

y = -y-af + jf, ( 2 ) 

to approximate the motion of stars orbiting in a plane about the galactic center, where 
the dot denotes the differentiation with respect to the time. Its solution is chaotic for 
some initial conditions, such as 

x(0) = 0.56,y(0) = 0,x(0) = 0,y(0) = 0, 

as mentioned by Sprott [16J. Without loss of generality, let us use this chaotic system 
to describe the basic ideas of the CNS and to show its validity. 

It is well-known [21 El El HOI ED EH EE1 HE] that chaotic dynamic systems have the 
sensitive dependence on initial conditions (SDIC), i.e. small difference of initial condi- 
tions lead to great difference of numerical simulations for large time, so that long-term 
prediction of chaos is impossible. Since the traditional numerical simulations contain 
the unavoidable truncation and round-off errors at each time-step, which are not neg- 
ligible, all traditional numerical simulations of chaos are mixed with these numerical 
noises and thus are not "clean". Because these numerical noises are generally much 
larger than the micro-level inherent uncertainty of physical variables, the propaga- 
tion of such kind of physical uncertainty of chaotic dynamic systems has never been 
studied in details. 
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In order to gain reliable chaotic solutions in a long enough interval of time, Liao [9] 
developed a numerical technique with extremely high precision, called here the "clean 
numerical simulation" (CNS). Using the computer algebra system Mathematica with 
the 400th-order Taylor expansion for continuous functions and data in accuracy of 
800-digit precision, Liao [9] gained, for the first time, the "clean" numerical results 
of chaotic solution of Lorenz equation in a long interval < t < 1000 LTU (Lorenz 
time unit) with negligible truncation and round-off error. The basic ideas of the 
CNS are simple and straightforward. Since the order of Taylor expansion is very 
high, the corresponding truncation error is rather small. Besides, expressing all data 
in the accuracy of large-number digit precision, the small enough round-off error is 
guaranteed. Thus, as long as the order of Taylor expansion is high enough and the 
digit-number of data is long enough, both of the truncation and round-off errors 
can be much smaller than the micro-level uncertainty so that the propagation of 
micro-level uncertainty can be simulated accurately in a long enough interval of time. 
Currently, Liao's "clean" chaotic solution [9] of Lorenz equation is confirmed by Wang 
et al [17J, who used parallel computation with the multiple precision (MP) library: 
they gained reliable chaotic solution up to 2500 LTU by means of the lOOOth-order 
Taylor expansion and data in the accuracy of 2100-digit precision, and their result 
agrees well with Liao's one [9] in < t < 1000 LTU. This kind of reliable "clean" 
chaotic solutions and especially the CNS provide us a way to accurately investigate 
the propagation of micro-level uncertainty of chaotic dynamic systems, which may 
give us some new explanations to the essence of SDIC and the "butterfly effect" from 
the physical and statistic points of view, as shown below. 

As mentioned above, the basic ideas of the CNS are easy to understand. The first 
step of the CNS is to give the high-order Taylor formulas for numerical simulation. 
Let ( 

x ni Dn) and (x n , y n ) denote the position and velocity at the time t n — nAt, where 
At is a constant for time-step. Write 



where ao = x n ,bo = y n , a\ = x n ,bi = y n are known. Substituting them into the 
original governing equations of the Henon and Heiles system [7] and equaling the like 
power of 5t = t — t n , we have the recursion formula 




x 



(t) = ^a n {t-t n ) n , y(t) = J2 b n (t-t n ) n , t n <t<t n + At 



n=0 n=0 



n 



a n + 2 a kb n -k 
k=0 



(3) 
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for n > 0. Then, we have the Mth-order Taylor approximation 
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and 

M-1 

x n+1 « ^(A; + l)a fe+1 (At) fc , (6) 

fc=0 
Af-1 

2/n+i « ^(fc + l)6 fe+1 (At) fc (7) 

fc=0 

at the time t n +i = (n+ I) At. Besides, all data must be expressed in a large enough 
number of digit precision. In this article, for the sake of simplicity, all data are 
expressed in the accuracy of 2M-digit precision (we use the computer algebra system 
Mathematica here), as long as M, the order of Taylor expansion, is large enough. In 
this way, one gains the high-order numerical simulations of x(t) and y(t) step by step, 
with small enough truncation and round-off errors. 

Let T c denote the maximum time for given M and At, up to which the chaotic 
solution is "clean", i.e. without observable influence of the round-off and truncation 
errors, and thus is reliable. The second step of the CNS is to determine a good 
approximation of T c , an important characteristic length-scale of time for the CNS. 
Without loss of generality, we use here the Mth-order Taylor formula with At = 1/10 
and the data in accuracy of (2M)-digit precision. It is found that 

T c ss 32 + 32M, 

and 

T c w21-231og 10 H, 

i.e. 

| a | < io°- 875 ^ - 044Tc 

where a is the standard deviation of physical fluctuation of position and velocity 
in initial conditions. For examples of gaining such kind of formulas, please refer to 
Liao [9]. Thus, in order to gain "clean" reliable chaotic solution of the Henon and 
Heiles system [7] in the interval < t < 2000, we should use at least the 70th-order 
Taylor expansion (with At = 1/10) and the data in accuracy of 140-digit precision. 

Note that the numerical noises of each simulation given by the CNS are negligible 
in the interval < t < 2000, so that we can investigate the propagation of the inherent 
micro- level uncertainty of the position and velocity at t — 0, which are much larger 
than the numerical noise, as shown below. 

3 The micro-level uncertainty 

The kinetic status of a star is determined by its position and velocity. Traditionally, 
it is believed that the kinetic status of a star is inherently exact and the uncertainty 
of position and velocity come from the imperfect measure equipments. However, as 
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pointed out below, this traditional idea is suspectable: the position and thus the 
velocity of a star are inherently uncertain. Thus, let us consider the initial conditions 

ar(0) = x + x Q , y(0) = y + y , x(0) =u + u , y(0) = v + v , 

where xo, yo, uo, vq are observable values of the initial position and velocity of a star, 
and xo,yo,uo,vo are the corresponding micro-level uncertain ones, respectively As- 
sume that (xo, yo, uo, vo) is exactly defined and the uncertain term (x , yo, uo, ^o) is in 
the normal distribution with zero mean and a micro-level deviation a = 10~ 60 . The 
reasons for the above assumptions are given below. 

It is a common belief of the scientific society that the microscopic phenomenon 
are essentially uncertain and random. To show this point, let us consider some typical 
length scales of microscopic phenomenon widely used in modern physics. For example, 
Bohr radius 

h 2 

r = « 5.2917720859(36) x KT 11 (m) 

m e e 1 

is the approximate size of a hydrogen atom, where h is a reduced Planck's constant, m e 
is the electron mass, and e is the elementary charge, respectively. Besides, Compton 
wavelength L c = h/(mc) is a quantum mechanical property of a particle, i.e. the 
wavelength of a photon whose energy is the same as the rest-mass energy of the 
particle, where m is the rest-mass of the particle and c is the speed of light. It is 
the length scale at which quantum field theory becomes important. The value for the 
Compton wavelength of the electron is 

L c » 2.4263102175(33) x 10~ 12 (m). 

In addition, the Planck length 



l P = J— ss 1.616252(81) x 10" 35 



c 3 



where c is the speed of light in a vacuum, G is the gravitational constant, and h is the 
reduced Planck constant, is the length scale at which quantum mechanics, gravity and 
relativity all interact very strongly. Especially, according to string theory [14], the 
Planck length is the order of magnitude of the oscillating strings that form elementary 
particles, and shorter length do not make physical senses. Besides, in some forms 
of quantum gravity, it becomes impossible to determine the difference between two 
locations less than one Planck length apart. Therefore, in the accuracy of the Planck 
length level, the position of a star is inherently uncertain, so is its velocity. Note 
that this kind of uncertainty is inherent and has nothing to do with the Heisenberg 
uncertainty principle [6]. 

On the other hand, according to de Broglie [2j , any body has the so-called wave- 
particle duality, and the length of the so-called de Broglie wave is given by 



A-^A-f-V, (9) 



mv V 
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where m is the rest mass, v denotes the velocity of the body, c is the speed of light, 
h is the Planck's constant, respectively. Note that, the de Broglie's wave of a body 
has non-zero amplitude, meaning that the position is uncertain: it could be almost 
anywhere along the wave packet. Thus, according to the de Broglie's wave-particle 
duality, the position of a body is inherent uncertain. 

Therefore, it is reasonable for us to assume that the micro-level inherent fluctu- 
ation of position of a star shorter than the Planck length l p is essentially uncertain 
and/or random. 

To gain the dimensionless Planck length l p , we use the dimeter of Milky Way 
Galaxy as the characteristic length, say, du ~ 10 5 (light year) « 9 x 10 20 (m). 
Obviously, l p /dj^ ~ 1.8 x 10 -56 is a rather small dimensionless number. As mentioned 
above, two (dimensionless) positions shorter than 10 -56 do not make physical senses. 
Thus, it is reasonable to assume the existence of the inherent uncertainty of the 
position and velocity of a star in the normal distribution with zero mean and the micro- 
level standard deviation 10~ 60 , which is added to the observable values (x , y , u , v ) 
of the initial conditions. 

Generally speaking, the number 10 -60 is indeed rather small. The propagation 
of such a tiny inherent fluctuation of physical variables with chaos has never been 
investigated, since the numerical noises (i.e. truncation and round-off errors) of tradi- 
tional numerical simulations of chaos are much larger than the uncertainty in general. 
However, by means of the CNS with the 70th-order Taylor expansion and the data in 
accuracy of 140-digit precision, such kind of propagation can be accurately studied, 
for the first time, because 10~ 60 is a hugh number compared to 1CT 140 : the number 
10~ 140 /10 -60 = 10 80 is much greater even than dju/lp, the ratio of the dimeter of 
Milky Way Galaxy to the Planck length! 

4 Statistic property of chaos 

Without loss of generality, let us consider the case of the observable values 

x = 14/25, y = 0, m = 0, v = 
of the initial conditions, corresponding to a chaotic motion [15]. Let 



denote the sample mean and unbiased estimate of standard deviation of x(t), respec- 
tively, where Xi(t) is the zth sample given by CNS for a tiny random term (x , y , u , v ) 
with the micro-level deviation a = 10~ 60 in the initial condition, and iV = 10 4 is the 
number of total samples gained by the CNS mentioned above. 

The standard deviations a x {t) and o~ y (t) of x(t),y(t) are as shown in Figs. [Hand 
El respectively. Note that there exists an interval < t < Td with w 1000, in 
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Figure 1: The standard deviation o x of x in case of x = 14/25, y = 0, u = 0, v = 
and the uncertain term (xo,yo,Uo,Vo) in the normal distribution with zero mean and 
a micro-level deviation a = 10~ 60 . 
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Figure 2: The standard deviation a y of y in case of xq = 14/25, yo = 0, uq = 0, vq = 
and the uncertain term (xo,yo,Uo,Vo) in the normal distribution with zero mean and 
a micro-level deviation a = 10~ 60 . 
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Figure 3: The standard deviation a u of x in case of x = 14/25, y = 0, u = 0, v = 
and the uncertain term (x , y , u , v ) in the normal distribution with zero mean and 
a micro-level deviation a = 10~ 60 . 





Figure 4: The CDF of x', compared to the normal distribution (dashed line) with zero 
mean and the standard deviation of x' at t — 2000. Solid line: CDF of x' at t — 1500; 
symbols: CDF of x' at t = 2000. 
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Figure 5: The CDF of y', compared to the normal distribution (dashed line) with zero 
mean and the standard deviation of y' at t — 2000. Solid line: CDF of y' at t — 1500; 
symbols: CDF of y' at t = 2000. 

which a x (t) and cr y (t) are in the level of 10~ 14 so that one can accurately predict the 
position (x, y) of a star, even if the corresponding motion is chaotic and the initial 
condition contains uncertainty. Similarly, the velocity of the star can be also precisely 
predicted in < t < T d , as shown in Fig. [3] for the standard deviation a u {t) of x(t). 
Thus, when < t < T^, the behavior of the system looks like "deterministic" and 
predictable, even from the statistic viewpoint. When t > Td, the standard deviations 
of the position and velocity begin to increase rapidly, and thus the system becomes 
random obviously: the position (x,y) and velocity (x,y) of the star are strongly 
dependent upon their tiny inherent random terms (xq, jjo, Uq, Vq) of the initial condition 
in the micro- level accuracy of 10~ 60 . In other words, due to the SDIC of chaos, 
the unobservable micro-level inherent uncertainty of the position and velocity of a 
star transfers into the macroscopic randomness of the motion. So, chaos is a bridge 
from the micro-level uncertainty to macroscopic randomness! Therefore, the micro- 
level inherent uncertainty of the position and velocity is an origin of the microscopic 
randomness of motion of stars in our universe. This provides us a new physical 
explanation for the SDIC of chaos. For this reason, each "big bang" [12J will crease a 
completely different universe! 

Besides, it is found that the standard deviations of the position and velocity be- 
come almost stationary when t > T s , where T s pa 1300, as shown in Figs.[T]to[3j Thus, 
when Td < t < T s , the system is in the transition process from the "deterministic" 
behavior to the stationary randomness. It is interesting that the stationary standard 
deviations of x(t) and y(t) are about 1/3, and their stationary means < x > and 
< y > are close to zero. It means that, due to SDIC of chaos and the micro- level 
inherent uncertainty of position and velocity, a star orbiting in a plane about the 
galactic center could be almost everywhere in the galaxy at a given time t > T s . 
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Write the fluctuations x' = x— < x > and y' = y— < y >. The stationary 
cumulative distribution functions (CDF) of x',y' are almost independent of time, as 
shown in Figs. H] and EJ Besides, the stationary CDF of the fluctuation x' is rather 
close to the normal distribution with zero mean and the standard deviation of x', as 
shown in Fig. HI But, the stationary CDF of the fluctuation y' is obviously different 
from the normal distribution, as shown in Fig. |5j 

Similarly, we investigate the influence of the observable values (xq, yo, Uq, vq) and 
the standard deviation a of the uncertain terms (x , y , u , v ) in the initial condition. 
It is found that T d decreases exponentially with respect to cr. Besides, the stationary 
means and standard deviations of x, y, x, y, and the CDFs of x' and y', are inde- 
pendent of the observable values (x , y , u , v ). Thus, when t > T s , all observable 
information of the initial condition are lost completely. In other words, when t > T s , 
the asymmetry of time breaks down so that the time has a one-way direction, i.e. the 
arrow of time. So, statistically, the Henon-Heiles system has two completely different 
dynamic behaviors before and after T^: it looks like "deterministic" and predictable 
without time's arrow when t < T^, but thereafter rapidly becomes obviously random 
with arrow of time. 

5 Conclusions and discussions 

In this paper, a fine numerical algorithm, namely the "clean numerical simulation" 
(CNS), is proposed to accurately simulate the propagation of micro-level uncertainty 
of chaotic dynamic systems. The chaotic Henon-Heiles system is used as an example 
to show the validity of the CNS. 

Based on the assumption that the micro-level inherent fluctuation of position 
of any a body is essentially uncertain and/or random in the accuracy of micro- level 
shorter than the Planck length l p , our CNS simulations indicate that, due to the SDIC 
of chaos, the micro-level inherent uncertainty of position and velocity of a star transfers 
into macroscopic random motion. Thus, the macroscopic randomness of motion of 
stars has a close relationship with the microscopic uncertainty of position and velocity: 
chaos is a bridge from the micro-level uncertainty to macroscopic randomness. This 
kind of relationship gives an explanation and origin of the random motion of stars in 
the universe. It also provides us a new physical explanation to the SDIC of chaos. 
Since chaos and the microscopic uncertainty of position and velocity of a star are 
unavailable, each "big bang" [12] will create a completely different universe. 

Note that the Planck length l p in the assumption is not the key point. Replacing 
l p by 10~ 30 x l p , we gain qualitatively the same conclusions. Besides, our assumption 
is consistent with the Heisenberg uncertainty principle [6] 

Ax A(mv) > h, 

where h is the Planck constant, which indicates that 

h h 
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say, it is impossible to measure a position more accurate than h/(mc), since the un- 
certainty of velocity Av is not greater than the light speed, according to the relativity 
theory |4J. So, according to the Heisenberg uncertainty principle [6J, it is impossible 
to measure the position of a star (with the mass m) more accurate than h/(mc). 
Thus, since m and c are finite, it is impossible to measure a position with an arbitrary 
degree of accuracy and certainty, forever. So, considering the string theory [TJ] and de 
Broglie's wave-particle duality, we have many reasons to make the assumption used in 
this article: the position of a star is inherently uncertain, and such kind of uncertainty 
has nothing to do with measurement done by human being. 

Strictly speaking, the micro-level uncertain fluctuation of the position and velocity 
of a star exists not only at t = 0, but also in the whole interval of time. Adding the 
micro-level fluctuation to the calculated results at each time step, we gain the same 
stationary means, standard deviations and CDFs of x,y, and the same conclusion 
that the inherent micro-level uncertainty of the position and velocity transfers into 
the macroscopic randomness of motion of a star. 

Note also that the Henon-Heiles system is a model for a star orbiting in a plane 
about the galactic center. Any a mathematical model is more or less a kind of approx- 
imation of physical phenomena. However, one can gain the qualitatively same conclu- 
sions, even if the Henon-Heiles system is replaced by other better models someday, as 
long as their solutions are chaotic in some cases. Due to the Newton's gravitational 
law, the motion of stars are governed by nonlinear dynamic equations so that chaotic 
motions of stars are unavailable [TJ. Thus, if the equations of the Henon-Heiles system 
are replaced by those for the famous three-body problem [13] . one can gain the qual- 
itatively same conclusions. So, our conclusions reported in this article have general 
meaning. 

Traditionally, the deterministic means that the evolution of solutions is fully de- 
termined by initial conditions without random or uncertain elements involved. Strictly 
speaking, such kind of deterministic models do not exist at all for physical phenom- 
ena, since no one can ensure that there exist no random and uncertain elements in 
initial/boundary conditions. So, strictly speaking, the Henon-Heiles system is not a 
deterministic model even from the traditional viewpoint. Similarly, the three-body 
problem is also not a deterministic one, since the inherent uncertainty of position 
and velocity in its initial condition is unavailable, as mentioned above. Thus, strictly 
speaking, the traditional "deterministic chaos" does not exist for models having phys- 
ical meanings: only a few pure mathematical models without physical meanings are 
"deterministic chaos". In this article, we illustrate that a nonlinear dynamic system 
with chaos can be predicted with a rather high-level accuracy in a long interval of 
time and thus behaves like "deterministic" and predictable, even from the statistic 
viewpoint. This might provide us a new definition of the "deterministic", which is 
more practical than the traditional ones. 

It is traditionally believed that the long-term prediction of chaos is impossible, 
mainly due to the impossibility of the perfect measurement of initial conditions with 
an arbitrary degree of accuracy. This is the traditional explanation to the SDIC of 
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chaos. Here, we provide a new explanation for the SDIC of chaos from the physical and 
statistic viewpoint: the observable initial conditions are assumed to be exact, but the 
inherent micro-level uncertainty of the physical variables are propagated into macro- 
scopic randomness. Different from the traditional explanation of the SDIC, which 
focuses on the measurement, the new explanation emphasizes the inherent micro-level 
uncertainty of physical variables and its propagation with chaos. Besides, it should 
be emphasized that such micro-level uncertainty of physical variables was completely 
inundated with the numerical "noises" of the traditional numerical methods, and thus 
have never been studied in details. This shows the validity and potential of the CNS 
to precisely simulate complex physical phenomena with the SDIC, such as weather 
prediction and turbulence. 

Finally, it should be emphasized that, the macroscopic randomness is traditionally 
believed to have no relationships with the micro-level uncertainty. However, our work 
illustrates that they have a close relationship: chaos is a bridge from the micro-level 
inherent uncertainty to macroscopic randomness! 
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